Optical Absorption Study by Ab initio Downfolding Approach: Application to GaAs 
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We examine whether essence and quantitative aspects of electronic excitation spectra are correctly captured by 
an effective low-energy model constructed from an ab initio downfolding scheme. A global electronic structure 
is first calculated by ab initio density-functional calculations with the generalized gradient approximation. With 
the help of constrained density functional theory, the low-energy effective Hamiltonian for bands near the Fermi 
level is constructed by the downfolding procedure in the basis of maximally localized Wannier functions. The 
excited states of this low-energy effective Hamiltonian ascribed to an extended Hubbard model are calculated 
by using a low-energy solver As the solver, we employ the Hartree-Fock approximation supplemented by the 
single-excitation configuration-interaction method considering electron-hole interactions. The present three- 
stage method is applied to GaAs, where eight bands are retained in the effective model after the downfolding. 
The resulting spectra well reproduce the experimental results, indicating that our downfolding scheme offers a 
satisfactory framework of the electronic structure calculation, particularly for the excitations and dynamics as 
well as for the ground state. 



I. INTRODUCTION 

First-principles electronic-structure calculations based on 
density-functional theory '(DFT) within the local density ap- 
proximation (LDA) or the generalized gradient approxima- 
tion (GGA) for the exchange-correlation (XC) functional have 
opened a way to predict ground-state properties of various 
materials without introducing ad hoc parameters. However, 
there exist serious problems in which the DFT fails even qual- 
itatively. Typical examples are found in strongly-coiTelated 
electron systems such as the genuine Mott insulator, where an 
insulating gap opens in partially filled bands solely owing to 
the strong local electron-electron repulsion.^ The DFT with 
LDA/GGA often predicts metals for these systems,^ indicat- 
ing the fact that the XC functionals based on LDA/GGA do 
not correctly capture the local coiTelation in real space. Other 
typical example is found in dynamics and excitation spec- 
tra of electrons, in which many-body correlation effects are 
also essential.'*'^ Even semiconductors, being supposed to be- 
long to weakly-correlated electron systems in the ground state, 
may have highly-degenerate excited states arising from the lo- 
cal electron coiTelation effects, and thus the single-particle 
approximations such as the Kohn-Sham^ and Hartree-Fock^ 
schemes break down in general. It is well known that incor- 
porating two-particle interactions between electrons and holes 
generated by the excitation is crucial in describing the elec- 
tronic structure at low-energy levels. A typical example is 
found in excitonic excitations.^'^'^"' ' 

To treat these excitations properly, we clearly need to go 
beyond the single-particle theory, while a full ab initio calcu- 
lation taking into account the many-body correlation effects 
is practically intractable. To go beyond the LDA/Hartree- 
Fock levels, we are required to develop a sufficiently accurate 
but efficient and practically feasible method. This challenge, 
so-called "beyond LDA/Hartree-Fock" problem has attracted 



growing interest. The GW method'-''^ has been devel- 
oped to incorporate self-energy effects basically on the level 
of the random phase approximation (RPA) while strong cor- 
relation and fluctuation effects beyond the RPA level require a 
more accurate and reliable treatment. Especially, an ab initio 
three-stage scheme has been rapidly developed by combin- 
ing two procedures, namely, LDA/Hartree-Fock framework 
and accurate low-energy solvers. '^"'^ The global electronic 
structure is first obtained by the LDA/Hartree-Fock scheme. 
In the next stage, one performs a bridging treatment, that is 
downfolding, ^^•^^'^'^ by eliminating the high-energy degrees of 
freedom leaving the low-energy effective model (Hamiltonian 
or Lagrangian) for local bases like Wannier functions. ^^'^^ 
The downfolding determines parameters for the effective low- 
energy model via first-principles calculations. The resulting 
low-energy model is, in the final stage, solved by low-energy 
reliable solvers such as dynamical mean field theory,'^ path- 
integral renormalization group, and/or various Monte- 
Carlo methods^^ developed for treating the coiTelation effects. 
Such a hierarchical three-stage scheme instead of a full ab 
initio calculation allows us to perform a first-principles and 
parameter-free prediction of the electronic structure of the 
strongly-correlated electron system within the present feasi- 
bility of computer. 

In this paper, we present theoretical studies on the ab initio 
downfolding scheme to assess the reliability for treating dy- 
namical properties. In our scheme, maximally localized Wan- 
nier functions are introduced as a basis function for represent- 
ing the model Hamiltonian. This basis offers computation- 
ally convenient choice, because this Wannier function can be 
computed with any basis functions (plane wave,^^ linearized 
muffin-tin orbital,^'' linearized augmented plane wave,^^ etc). 
Transfer parameters are evaluated by calculating Kohn-Sham 
matrix elements in this basis, and onsite/offsite interaction pa- 
rameters including screening effects are determined via con- 
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strained calculations. ' • 

In the three-stage scheme, the rehabihty of the downfold- 
ing procedure and the accuracy of the resulting model param- 
eters are crucially important. In particular, the reliability in 
describing dynamics and excitation spectra has to be critically 
tested. For this purpose, we make a critical comparison be- 
tween experimental results and computational results for the 
generated model. In the present study, we focus on optical- 
absorption properties. It is widely recognized and accepted 
in the literature that the optical absorption of solids, in partic- 
ular for semiconductors or insulators, is deeply affected by 
excitonic effects.'*'^ This effect originates from an effective 
Coulomb interaction between electrons and holes and there- 
fore is sensitive to the magnitude of the interaction parame- 
ters in the model Hamiltonian. To examine this effect through 
the present formahsm, we choose GaAs as a representative 
material exhibiting spectral enhancement due to the excitonic 
effect and calculate its optical spectra by taking account of the 
electron-hole interaction. There exist many experimental^ ''^^ 
and highly-accurate ab initio^^'^^ spectral data for this mate- 
rial. Therefore, our downfolding formalism and determined 
parameters can be critically tested by examining whether our 
model spectrum reproduces those data satisfactorily. 

This paper is organized as follows: In Sec. II we describe 
our downfolding procedure; we introduce a complete-neglect- 
differential-overlap model which is used as our target model 
Hamiltonian and describe computational details for determin- 
ing the model parameters. In Sec. Ill, to take into account 
the electron-hole interaction, we introduce a single-excitation- 
configuration-interaction framework for calculating an optical 
absorption. Efficient techniques to evaluate one-body velocity 
matrix elements needed in the spectral calculation, based on 
the Wannier interpolation scheme, is described in appendix. 
The calculated optical spectra are compared with the experi- 
mental results. Concluding remarks are given in Sec. IV. 



II. DOWNFOLDING PROCEDURE 
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FIG. 1: (Top) A global ab initio band structure of GaAs at [—15 
eV:+30 eV]. (Bottom) ab initio original (solid line) and interpolated 
(dots) bands. Energy zero is set to the top of the valence bands. The 
energy window is set to [—15 eV:-|-10 eV]. 



A. Global electronic structure by DFT 

The first procedure derives the global electronic band struc- 
ture by a conventional DFT scheme. The present scheme is 
based on ab initio density functional calculations with Tokyo 
Ab initio Program Package^^ developed by the condensed- 
matter-theory group in the University of Tokyo. With this 
code, band calculations have been performed within the 
generalized gradient approximation^^ to the exchange cor- 
relation functional, using a plane-wave basis set and the 
TrouUier-Martins norm-conserving pseudopotentials''^ in the 
Kleinman-By lander representation."*" The energy cutoff is set 
to 25 Ry, and a 15x15x15 fc-point samphng is employed to 
represent electronic structures of the system. The resulting 
global band structure of GaAs at an energy region [—15 eV:+ 
30 eV] is illustrated in the top panel of Fig. 1 



B. Complete neglect differential overlap model 

Now we go onto the second stage and start the derivation 
of an effective low-energy Hamiltonian by the downfolding 
procedure. Before going to the downfolding itself, we first 
specify the form of the final effective low-energy Hamilto- 
nian. In the present downfolding procedure, we will make 
several approximations by simplifying the low-energy effec- 
tive Hamiltonian, in which we expect that the approximations 
do not alter the optical spectra seriously. The first approxima- 
tion is to consider only the diagonal Coulomb interaction and 
ignore the offdiagonal part (exchange interaction) in the low- 
energy effective model, which results in an extended Hubbard 
model or, in other words, the complete-neglect-differential- 
overlap (CNDO) model. The CNDO Hamiltonian has origi- 
nally been introduced by Pople^^ to study electronic structures 
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of small organic molecules and, to date, has been extended to 
study various properties of complicated systems ranging from 
transition-metal compounds^'* to proteins^^ and DNA.^^ A re- 
markable property of this Hamiltonian is that it considers all 
the degrees of the freedom of valence electrons of the system, 
which allows describing the individual characters of the real 
material. 

The crystal CNDO Hamiltonian H consists of a one-body 
part Ht and an interaction part He'- 

n = Ht + nc. (1) 

The one-body part Tit is given by 

+ EEEwj(^'-^)<fl<,fl'}. (2) 

RR' ij iiv 

where a^J^ (^^ir) is a creation (annihilation) operator of a 
valence electron with spin a in /U-type localized basis centered 
at ith site in lattice R. As mentioned in the introduction, we 
use the maximally localized Wannier function (MLWF) as the 
basis function for representing the CNDO Hamiltonian. In the 
present study of GaAs, there are eight WF's in the primitive 
cell, where the first four belong to a Ga site, and the remain- 
ing four belong to an As site. Thus, the index /i specifies four 
types of lobe directions (band indices) of the MLWF's, and the 
suffix i specifies the Ga or As sites, /^i and t^i^j {R—R') are 
the ionization potential and the transfer integral, respectively. 
Notice that the translational symmetry in the crystal is exphc- 
itly considered for matrix elements; /^jjt = I^i and t^iR^jR' 
= t^iiyj {R' — R) for any R and R'. 
The interaction part He is written as 

'^c = E{EEf^^^Mifl^iifl+EE^i^M'«^-«} 
+ EE^«^(^-^')(^iii-^i)(^.fl'-^.)- (3) 

RR' ij 

Here, iV^^ = a^lfla-,^, N^^r = E.^^iR' and NiR = 
^fiiR are the number operators, and Zi is the core charge. 
Ui and Ul are the onsite intraorbital and interorbital Coulomb 
repulsions, respectively. Vij{R — R') in the third term is an 



interatomic Coulomb repulsion, and it is assumed to be inde- 
pendent of the lobe directions /j. and i^. 



C. Parameterization 

We now describe the downfolding procedure and parame- 
terization for the CNDO model of Eq. (1). The downfolding 
consists of two parts. The first is the derivation of the kinetic 
energy part, where a tight-binding Hamiltonian Hi is derived. 
The second part is the derivation of the interaction part He- 
1. Kinetic Energy 

The tight-binding Hamiltonian Ht given in Eq. (2) is de- 
rived from the global band structure after eliminating higher- 
energy bands. This downfolding may be performed by the 
perturbation scheme. ^^ '^ The resultant band structure is nor- 
mally very close to the low-energy part of the original band 
structure and the difference is not discernible when the low- 
energy retained part is isolated from the eliminated high- 
energy bands. This means that the self-energy of the re- 
tained bands caused by the higher-energy eliminated electrons 
is negligible. Since such self-energy effects are smaller even 
in semiconductor systems, in this paper, we employ the low- 
energy part of the bands as the retained bands after the elimi- 
nation of the higher-energy bands. 

Now we retain eight bands near the Fermi level and con- 
struct Wannier orbitals from the retained band structure. To 
this end, ab initio MLWF's are constructed with the Souza- 
Marzari-Vanderbilt algorithm.^^ We set an energy window in 
the interval [—15 eV:+10 eV], which includes four valence 
and four conduction bands of the system. The resulting Ga 
and As MLWF's are displayed in the top and bottom panels 
of Fig. 2, respectively. We see that the Wannier functions are 
almost locaUzed at a single site and have an anisotropic char- 
acter due to an sp^ hybridization. To show the accuracy of 
low-energy band structures represented by the resultant WF's, 
we compare in the bottom panel of Fig. 1 original bands (solid 
line) with interpolated bands (dots) obtained by diagonalizing 
fc-space Kohn-Sham (KS) Hamiltonian matrices represented 
by the WF's. We see a good agreement between the original 
and interpolated bands in the energy window. Ionization po- 
tential Ifi^i and transfer integral (i?) are extracted from 
the matrix elements of the one-body KS Hamiltonian /iks in 
the basis of the MLWF's \wniR)as 



Ifii = {w^,io\hKs\w^io) and tp,://jf-R) = ("',-'o|^'ks|w',-,,k) 



(4) 



respectively. Here R = lai + ma2 + nas with —7 < 2. Interaction Energy 

{I, m, n) < +7, and {ai , 02 , as} are primitive lattice vectors. 

We next derive the interaction parameters Ui, 17/, and 

Vij{R) for the low-energy model. In the original CNDO 
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FIG. 2: (Color online) Maximally localized Wannier functions of Ga 
(top) and As (bottom). The amplitudes of the contour surface are 
+0.5/v'u (blue) and (red), where v is the volume of the 

primitive cell. The shaded sheet represents a 3 x 3 x 3 fee lattice and 
Ga and As nuclei are illustrated as gray and blue dots. 



framework, the Coulomb interaction between electrons has 
long-ranged tail as scaled by 1 /r in the limit r ^ oo with 
r being a distance between two electrons. We are, however, 
interested in the electronic interaction in condensed phase for 
which the Coulomb interaction is effectively screened. In fact, 
the dielectric constant of GaAs is rather high; eo = 10.6 
experimentally,'*^ indicating that the Coulomb interaction de- 
cays as 1/ (ear) in the limit r ^ oo. We thus employ an 
approximation for this interaction; we keep only onsite inter- 
actions Ui and C/j' and the nearest-neighbor interaction V. 

These parameters are determined with a constrained DFT 
framework following a "hopping-cutoff ' treatment."*^ As the 
basic strategy, one first kinetically decouple a specific site 
from the rest of the system, thus leaving this site isolated 
as the so-called atomic limit. This decoupling treatment is 
achieved by switching off the hoppings between the Wannier 
orbitals at the specific site and all the other orbitals, where we 



identify the hoppings with the off-diagonal Kohn-Sham ma- 
trix elements in the representation of the Wannier functions. 
With such a hopping cutoff, the standard constrained total en- 
ergy calculations are performed to generate a potential energy 
surface with respect to constrained parameters such as occu- 
pancies of the Wannier orbitals belonging to the decoupled 
site. The interaction parameters obtained with quadratic fit- 
ting to the resulting potential energy data include screening 
effects ascribed to the relaxation of the valence electron den- 
sity around the decoupled site. In the present case, the proce- 
dure for determining the interaction parameters is somewhat 
complicated because of a large number of the interaction pa- 
rameters to be derived. So, in the practical work, we divide 
the treatment into two steps; the determination of U and U' 
and the subsequent step for determining an offsite parameter 
V. 

The basic strategy for obtaining U and U' is to generate 
potential-energy-data sets with respect to two types of the con- 
strained parameters; (I) the first potential-energy data are ob- 
tained from the constrained calculations with respect to occu- 
pancy Qfj^j of a specific Wannier orbital /i at the site /, and (II) 
another data are obtained from the constrained calculations for 
a site occupancy Qj defined as the total amount of the orbital 
occupancies belonging to this site; g^/. The curvature of 
the first potential energy curve plotted as a function of the or- 
bital occupation (7^/ gives an estimate of the onsite intraorbital 
interaction for the orbital fi at the site /, while the curvature 
of the second data represents the averaged value over the on- 
site intraorbital/interorbital interactions. From quadratic fit- 
ting to the mixed two data, we can determine the U and U' 
parameters reasonably (see below). Though the types of the 
constraints are different in the two calculations, the method 
itself is the same, so, here, we describe only details for the 
constrained calculation with respect to the single orbital occu- 
pancy. 

The practical calculation proceeds as follows: We first con- 
sider a 3 X 3 X 3 fee supercell containing 54 atoms,**' and 
choose the central Ga site placed at the origin as the decou- 
pled site. (Here, we describe only the Ga case, but a parallel 
treatment can be applied to the U and U' determination for the 
As site.) We next introduce a cutting operator Acut to switch 
off the hopping integrals connecting the four Wannier orbitals 
of this site to the other Wannier orbitals, 

Acut = -Pioh-KsPw - PwhKsPio + PiofiKsPio 
+ X! I^Vo)^M/(^M-fo| , (5) 

where Iiks is a one-body KS Hamiltonian, and Pw is a pro- 
jector onto the total Wannier orbitals, 

X i fi 

Here, X is a lattice vector denoting a supercell, the suffix i 
specifies the sites in the supercell, and fi stands for the band 
index of the Wannier orbital. P/o in Eq. (5) is a projector onto 
the Wannier orbitals belonging to the decoupled / site in the 
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home cell (X = 0), 

-Pro = ^ \w^tio){w^io\ . 



(7) 



With the cutting operator in Eq. (5), we define the constrained 
total energy as 



It^^ fak\{Wij.Io\(l)ak)\'^ - q^IO 



N 



(8) 



r 



Here, J^[p\ is a usual density functional with a total charge 
density p{r) = ^ Y.ak fak\(l)ak{r)\^ with N being the to- 
tal number of k points, the first term in the bracket [• • • ] is 
the definition itself for the orbital occupancy of the discon- 
nected Wannier orbital |w;^/o), and A^/o is the Lagrange mul- 
tiplier associated with the constraint to fix the orbital occu- 
pancy at q^io. A functional derivative of .Ectot with respect to 
the Bloch orbital (j)ak leads to the following constrained KS 
equation, 



hxs + AcMt + A^/o|'u;^/o)('U'^/o| 



|0ak) = eak\4>ak), (9) 



where the third term in the left hand side is an additional po- 
tential due to the occupancy constraint. For numerical details 
for solving the equation, readers are referred to Ref. 29. By 
solving the equation, we generate constrained potential en- 
ergy data (we refer to these data as DATA I), and plot them as 
functions of the orbital occupancy q^io and the site occupancy 
Qio = J2fj, Quia- In parallel to this treatment, the constrained 
calculations for the site occupancy are performed, where there 
is a small modification in the constrained KS equation (9); 
the additional potential is changed to A/o J^fi \'^fiio){'Wnio\. 
We again monitor the calculated constrained total energies 
as functions of g^/o and Qjo (DATA II). With the resulting 
DATA I and II, we perform quadratic fitting of the following 
function around g^/o = q^io and Qio = Qiq: 

1 - 2 

/ {Qio, QixIo) = 2^Ga. {Qio - Qio) 

+ 2 (C^Ga - f^Ga) {QiO " Qio) {QhIO - 



-I- 2 (f/ca - UqJ {q^io - q^io) 



(10) 



where g^/o and Qiq are equilibrium occupancies taken from 
the global band structure with no additional potential (A^/o = 
0). The form of the fitting function in Eq. (10) is derived 
by exploiting the character of the fourfold degeneracy of the 
Wannier orbitals {w^/o} (Appendix A). The Uoa and Uq^ 
values thus obtained are 2.39 eV and 2.17 eV, respectively, 
which are largely reduced from the bare interaction values 
t/g^ = 9.25 eV and U^^ = 7.89 eV. The same procedure 
is applied to the U and U' determinations of the As site. It 
was found to be Uas = 2.71 eV and J/^^ = 2.09 eV. The cor- 
responding bare values U^^ and U'^^ axe 11.45 eV and 9.80 
eV, respectively. 



We next describe the determination of the offsite parame- 
ter V. The interaction depends on the relative configuration 
between the Wannier orbitals. For example, let us consider a 
configuration formed by a pair of the Ga Wannier orbital and 
the As Wannier orbital, where these orbitals face along the 
covalent bond of the two atoms (we call it facing configura- 
tion). An example for the facing configuration can be found in 
the two Wannier orbitals displayed in the top and bottom pan- 
els of Fig. 2. Obviously, the strength of the Coulomb repul- 
sion in the facing configuration is relatively large compared to 
that in the other configurations. It should be noted here that 
the V parameter affects renormalized transfer integrals [see 
Eq. (12) in Sec. II D for the explicit form of the renormalized 
transfer integral]. Since bonding and anti-bonding orbitals are 
formed by the Wannier orbitals in the facing configuration, 
the bonding and anti-bonding spUtting, i.e., band gap itself, is 
dominated by the magnitude of a renormalized transfer inte- 
gral between the Wannier orbitals in the facing configuration. 
Thus, the V value in this configuration is crucial for an accu- 
rate description of the low-energy band structure and optical 
excitations. The contributions from the intersite interaction in 
the other configurations are much smaller in optical response. 
The result would not change when we shghtly overestimate 
them by V in the facing configuration, because the renormal- 
ized part to the transfer integral appears as the product of V 
and a density matrix [see Eq. (12)] and it was found that the in- 
tersite density-matrix elements are almost zero except for that 
of the facing configuration in the present GaAs case. There- 
fore, we calculate the V value in the facing configuration and 
employ it as the V value of the CNDO model. 

The actual determination of the V parameter proceeds as 
follows: We first choose two decoupled sites (the Ga site 
placed at the origin and the neighboring As site being in the 
[111] direction). The similar cutting treatment to Eq. (5) 
but extending the single-site-cutting formalism to the double- 
sites-cutting formalism is applied for this purpose. Then, we 
perform constrained calculations by imposing a constraint that 
two occupancies of the Wannier orbitals \w^io) and \wujo) 
in the facing configuration are kept at q^iQ and q^jo, re- 
spectively. We then draw the two-dimensional potential en- 
ergy surface with respect to the constrained parameters q^io 
and q^jQ, and perform a fitting of the quadratic function 

5C^Ga(a/:i/0 - QiJilof + ^UAs{qvJO - Qvjof + ^(O'/ii/O - 

qtiio){<li^Jo - Qujo) to the potential energy data. This fitting 
is performed by fixing [/ca and Uas at predetermined values 
in the preceding U and U' determination; we treat only V 
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as a single fitting parameter to avoid fitting errors and uncer- 
tainties. Tiie V value thus determined is 0.71 eV, where we 
again see the large reduction from the bare interaction value 
V° = 7.35 eV. The resulting interaction parameters are sum- 
marized in TABLE I. 

TABLE I: Interaction parameters determined in the present down- 
folding procedure. The energy unit is eV. 







Uas 




V 


2.39 


2.17 


2.71 


2.09 


0.71 



In general, the downfolded model contains an energy de- 
pendence in the interaction because the screening by the high- 
energy electrons necessarily causes a non-Markoffian and re- 
tardation effect. Such an energy dependence in the screened 
Coulomb interaction W{u}) is not described by the effective 
Hamiltonian. However, in the low-energy region, the screened 
Coulomb interaction is normally saturated to a constant and 
represented by the hmiting value at a; = 0. The constrained 
scheme is roughly regarded as the procedure to obtain this 
u! = limit. The effect of the larger (less screened ) interac- 
tion W at larger u as well as the effective interaction arising 
from virtual transition to eliminated bands can be accounted 



by the further consideration of the self-energy effect on the 
low-energy part.'^ For a wide band system such as GaAs, 
however, this effect may not be large^" and we thus ignore 
this effect. 



D. Hartree-Fock approximation 

In the electronic structure calculation of a typical semicon- 
ductor, GaAs, we expect that the strong correlation effect does 
not appear in the ground state and the Hartree-Fock (HF) ap- 
proximation provides us a reasonable result, although the ex- 
citation spectra are reliably determined only through the ac- 
count of the correlation effect more accurately. In this sec- 
tion, we consider how the HF solution for the ground state is 
calculated. 

The CNDO HF Hamiltonian in the basis representation of 
MLWF's is written as 

cr RR' ij IJ.V 

where -F^«;,yj (-R' — R) is the Fock matrix or the renormalized 
transfer matrix by the interaction part Tic of Eq. (3) and the 
matrix element is given by 



iQ^iO)-Zi)--{Q^i^i{0)-2) Ui 

+ [Qk{R" -R)- Zk] Vik{R" - R) 

R"k 

Wj(o)-^QMi.,(o)Vi,(o) 

tij,ivj{R' — R) — -Qfj,ivj{R' — R)Vij{R' — R) 



{R = R',i = j,iJ, = I'), 

(R = R',i = j, ji ^ v), 
(R = R,i^j), 
(R + R'), 



(12) 



Here Q{R — R) is the density matrix and the matrix element 
is given by Q^,i^j{R - R) = ($HF|Ea «^lR<jR' |*hf), 
where |<i>HF) is the HF ground state. The HF Hamilto- 
nian can be diagonaUzed with the Bloch orbital, Z*^^ = 

(l/V^v) E^i C-,{k) T.Re"'-"-0'l\w where a and fe are a 
band index and a wave vector, respectively, and N is the to- 
tal number of the unit cells in the system. The coefficients 
{C'^i{k)} are determined by solving the following Hartree- 
Fock equation, 



with 



R 



(14) 



With the resulting C^- (fe), the real-space density matrix is cal- 
culated by 



N 



R 



N 



(15) 



with 



Y F,i^j{k)C^j{k) = e„feC« (fe) (13) 



Q^,.,(fe) = 2^C«(fe)C«;(fe). 



(16) 
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by 



The CNDO total energy with the HF approximation is given where 



Etot = («'hf|W|$hf) 

~ 9 ^ ^ ^ QfJ,wj (fc) 



R 



(18) 



(17) The matrix element of the core matrix H"!^]^, (i?) is written as 



TTCOrO / TJN 



I^i-{Zi-l)Ui-^ZjVij{R') 
R'j 

^fxivj (-^) 



iR = 0,i = j,lj. = v), 
(otherwise). 



(19) 



r 



The actual CNDO HF calculation proceeds along the 
schematic diagram shown in Fig. 3. The initial density ma- 
trix is made from diagonalizing the core matrix H'^°'^'^{k) of 
Eq. (18). Since the Fock matrix F{R} in Eq. (12) depends on 
the density matrix Q{R) of Eq. (15), the HF equation (13) is 
solved self-consistently with an iterative procedure. To check 
the convergence, we monitor a density-matrix difference, 



EE 

ij JUL/ 



(0) 



(20) 



and a total-energy difference. 



(i-1) 



(21) 



where an upper suffix i stands for the number of the iteration 
step. The self-consistency condition we employ is satisfied 
when < IQ-^ and SE^^^ < 10-^ (a.u.). 

We show in Fig. 4 the self-consistent CNDO band struc- 
tures (solid fine), together with the ab initio interpolated band 
(dots). We see a rigid band shift in the conduction band. This 
rigid band shift is attributed to the renormalization of the inter- 
action part He in Eq. (3) into the one-body part; the so-called 
self-energy correction considered within the HF framework 
[see Eq. (12)]. We note that the trend of the rigid band shift is 
basically the same as the quasiparticle band shift observed in 
the GW calculation for semiconductor. ^^''^ 



III. ELECTRONIC EXCITATION 

In the previous section, we have described the procedure for 
the downfolding. We now start calculating physical quantities 
using the downfolded model. The purpose of this sections is to 
examine the rehabihty of the model obtained by the downfold- 
ing. In particular, we highUght whether the model gives us a 
reliable excitation spectra and dynamical properties. For this 
purpose, we calculate optical-absorption spectra, based on a 



configuration-interaction (CI) treatment considering electron- 
hole interactions, and compare the computational result with 
the experiments. 



A. Single excitation configuration interaction 

For electrons in sohds, the number of configurations gen- 
erated by the electronic excitations is in principle infinitely 
large, while the configurations capable of practical compu- 
tations are limited. Since we are interested in optical pro- 
cesses, we consider here only the single-excitation (SE) con- 
figurations which play a primarily important role in the hnear- 
absorption process, because the SE configurations directly 
couple with the HF ground state via an electric dipole opera- 
tor. The calculation at the SECI level in fact takes into account 
electron-hole interactions; the so-called excitonic effect in the 
spectrum. 

The SECI many-body wave function with a wave vector K 



IS written as' 



44 



u I/I 

fear 

where | ^ ^^k^^) is a spin-singlet SE configuration given by 

I'Kf'') = -^idlUKdl, + dtU^di,)\<^nF). (23) 

Here (dak) ^ creation (annihilation) operator of the 
Bloch electron in a virtual r (occupied a) band with spin a 
and a wave vector k. The CI coefficients {(7^^} in Eq. (22) 
and the excitation energy AE^k are obtained by solving the 
following CI equation, 

occ vir 

E E E ^^^rk.bsk' ^b^' = ^EgK C^fk (24) 
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Diagonalization ofH'^'^ky, Eq. (18) 



Construction of F(/f); Eq. (12) 



i 



Fourier transformation to F(ky, Eq. (14) 



1 



Diagonalization of F(A;); Eq. (13) 



i 



Construction of Qik); Eq. (16) 



i 



Fourier transformation to QiR); Eq. (15) 



i 



Total energy calculation; Eq. (17) 



i 



Convergence check; Eq. (20) and (21) 



No 



Yes 



Convergence 



FIG. 3: Schematic diagram for CNDO-HF band calculation 



with 



,K 



= {Crk+K - eakjSkk'SabSrs 

+ 2(rk+Kak\bk'Sk'+K) 
- {rk+KSk'+K\bk'ak), 



(25) 



where £Jhf is the HF ground-state eigenenergy for the many- 
body HF Hamiltonian Hhf in Eq. (11); TYhfI^hf) = 
-EhfI^hf). For the CNDO model, the second term in 
Eq. (25), called the exchange term, is calculated as 

{rk+Kak\bk'Sk'+K) = ^^C;:(fe + K)C;,(fe)C^*(fc') 

IJ,i vj 

X ClAk' + K)V^i.j{K), (26) 
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X W 



FIG. 4: Self -consistent GaAs CNDO-HF band structure (solid line) 
and ah initio interpolated band structure (dots). Energy zero is set to 
the top of the valence bands. 



evaluated by 

{I'k+K Sk'+K I bk' ak 



X C^*(fe')C« (fe)V^i.j(fe - fe'),(27) 



with 



and 



R 



(28) 



Ui {R = 0,i = j, ^ = ,y), 

(R = 0, ?: = J, ^ ^ v), (29) 



Vij {R) (otherwise). 

Eqs. (26) and (27) imply the repulsive-exchange and 
attractive-Coulomb interactions between an electron in the r 
and s bands and a positive hole in the a and b bands, with the 
total wave vector being kept constant [(A; -I- K) — k = K]. 

The structure of the present CNDO-HF-SECI equation (24) 
is basically the same as that of the Bethe-Salpeter equation 
for two-particle Green's functions. ^'^"'^^ A difference is that 
the former electron-electron interaction in the two-electron in- 
tegrals of Eqs. (26) and (27) is represented by Vfj,wj {R) de- 
termined via the constrained scheme (see Sec. II C), while 
the latter interaction is represented by the screened Coulomb 
interaction evaluated with the random phase approximation. 
Computationally, we note that four-center Coulomb integrals 
in the original exchange/direct terms are reduced to two-center 
Coulomb integrals including only the sites i and j because of 
the CNDO approximation.^^ Therefore the computational cost 
for the matrix evaluation is much smaller in our CI calcula- 
tion. The most time consuming step is the diagonalization of 
the CI matrix , which is scaled as the third power of the 



and the last term in Eq. (25), refereed to as the direct term, is dimension of the A^ ; (N ■ Nocc ■ ^vir 
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B. Optical absorption 

We next describe the expression for the SECI optical ab- 
sorption, which is given as the imaginary part of the macro- 
scopic transverse dielectric function,^ 



e2(w) = A/" 



(^'eJt=o|X|$HF) 



<5(A^eif=o-w),(30) 



where X. — Vi is a many-body position operator of elec- 
trons, and the normalization constant M is determined via the 
sum rule"*^ 



(31) 



with LUp being the plasma frequency of the system. Substitut- 
ing Eq. (22) into Eq. (30) and noting the commutation relation 

[Hhf,X] = Y.i ^HF(i), ri = - Y,i d/dr,, we obtain 



e k a 



eK- 

^ ark 



-_0 {(t>rk\d / dr\(t)ak) 



^rk ^ak 



(32) 



where {(pak) is the Bloch state being the eigenstate of the Fock 
operator /ihf and the matrix element of d/dr can be calcu- 
lated with an interpolation scheme^^ based on the MLWF's 
from first principles (see appendix B). Theoretically, £2(0;) 
contains an electron-hole-interaction effect due to the pres- 
ence of the CI coefficients in Eq. (32). To see the electron- 
hole-interaction effect on the spectrum, it is convenient to 
compare with the spectrum obtained by the independent- 
particle approximation**^ (IPA) where the CI coefficients are 
neglected in Eq. (32) and the excitations are described just 
with optical transitions between independent hole and elec- 
tron states. 



occ vir 



H = -^EEE 



k 

X 6{erk 



{(l>rk\d / dr\(j)ak) 



^rk ^ak 



^ak 



(33) 



In the spectral calculations, we have chosen a fc-grid dif- 
ferent from the Monkhorst-Pack fc-grid used in the band 
calculations.'** The present fc-grid is generated as follows: We 
first make a uniform fc-grid in an 11 x 11 x 11 Monkhorst- 
Pack mesh and then slightly shift uniformly the sampling k in 
the direction of -O.OI61 - O.O262 + O.O363 with {61, 62, ^3} 
being basic reciprocal lattice vectors. The resulting fc-points 
are different from the high-symmetry directions of the crys- 
tal and therefore are not connected via rotational operations 
of the crystal with each other. This leads to a finer sampling 
for the spectral calculation. An unshifted grid corresponds 
to only 56 crystallographically different points, which are too 
few to achieve a good spectral resolution. On the other hand, 
the shifted grid leads to a grid of 1331 crystallographically 
different fc-points, which gives a good spectral resolution. 



We show in Fig. 5 the calculated SECI (thick red line) and 
IPA (thin green line) spectra. The closed blue and open blue 
circles denote experimental results. ^''^^ We see a clear con- 
trast between the SECI and IPA spectra; by considering the 
electron-hole interaction with the SECI method, the spectral 
intensity in the low-energy region (< 5 eV) is enhanced, thus 
reproducing the experimental results perfectly. The agree- 
ment is indeed somewhat surprising, when we consider the 
several simplified treatments employed here such as the re- 
duction of the electron-electron interaction to the extended 
Hubbard form. However, we emphasize that the nature such 
as the spectral enhancement observed in proceeding from IPA 
to SECI is consistent with highly-accurate full ab initio re- 
sults'"" obtained from solving the Bethe-Salpeter equations. 
In the context of the downfolding, these results strongly sup- 
port that our model construction by the downfolding described 
in Sec. II C offers a reliable description not only of the ground- 
state band structure of an insulator but also of the excitation 
spectra. 








2 4 6 8 

Energy (eV) 
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FIG. 5: (Color online) Calculated SECI (thick red line) and IPA (thin 
green line) optical absorption spectrum of GaAs. The closed blue 
and open blue circles represent experimental data taken from Ref . 3 1 
and Ref. 32, respectively. 

For the completeness, we perform a more critical and elab- 
orate assessment of the reliability of the downfolding; we 
examine a sensitivity of the spectra to choices of Hamilto- 
nian parameters. In the present analysis, we focus on the 
check of the reliability of the offsite interaction parameter 
V. In fact, the magnitude of this parameter is expected to 
crucially control the strength of the electron-hole interaction 
[see Eqs. (26) and (27)] and thus directly affect the profile of 
the optical spectra. We may calculate the excitation spectra 
by using choices of interaction parameters different from the 
downfolded realistic values. Thus, the reliability of the down- 
folding can be assessed by examining whether the spectrum 
obtained from the present downfolded Hamiltonian gives the 
best agreement with the experiment among wider alternative 



choices of the interaction parameters. For this check, we intro- 
duce a scaHng factor x to scale V to xV . With this definition, 
x = 1 corresponds to the original ab initio V value, while 
in the region a; > 1, the nearest neighbor electron-electron 
repulsion is artificially overestimated. In the practical calcu- 
lation, we monitor the values of x at 0.8, 1.0, and 1.2 and then 
perform the SECI calculations to obtain the optical spectra 
for each V value. (In the calculations, the onsite parameters 
are fixed at the ab initio determined values displayed in TA- 
BLE I.) 

We show in Fig. 6 the resulting dependence of the SECI 
spectra (red lines) on the scaling factor x. The blue circles 
denote measured data. We see a notable change in the spectra 
due to the parameter increase [(a) (b) (c)]; increasing 
X makes a blue shift and an intensity decrease in the calcu- 
lated spectra. We see that the spectrum at the downfolded 
choice (i.e., the case with x = 1.0) exhibits the best agree- 
ment with the experiments among all the choices. The down- 
folded value offers the most realistic and accurate choices as 
the model parameters, and thus we conclude that the optical 
excitation spectrum is correctly captured by the downfolded 
Hamiltonian. 

Finally, we remark a recent development for ab initio evalu- 
ations for the offsite V parameter Indeed, applications of the 
constrained schemes to the determinations of the offsite pa- 
rameter is quite limited in the literature compared with those 
for the onsite parameters. Recently, Aryasetiawan et a/'^ have 
proposed an RPA approach for calculating the interaction pa- 
rameters. They first calculate a real-space screened Coulomb 
interaction C/(r, r') by excluding the polarization formed in 
target bands contained in the model Hamiltonian and then 
evaluate the matrix element of U in the localized basis such 
as linearized muffin-tin orbitals and/or maximally-localized 
Wannier functions. This approach is able to derive all the off- 
site parameters. So, comparisons between the values obtained 
by the present constrained scheme and the values based on 
the RPA approach would contribute to making deeper under- 
standing for the distant interaction parameters. 



IV. CONCLUSIONS 

We have examined whether the effective low-energy Hamil- 
tonian derived from the downfolding procedure is able to de- 
scribe dynamics and excitation spectra in a proper way. The 
calculation is performed in the three-stage scheme. In the 
first stage, we calculate the global electronic structure from 
the density functional theory supplemented by the generalized 
gradient approximation. The high-energy degrees of the free- 
dom in the global electronic band structure are, in the second 
stage, eliminated by the downfolding scheme, which leaves 
only the low-energy bands near the Fermi level. In the present 
example of GaAs, we retain up to 25 Ry for the calculation 
of the global electronic bands, while the downfolded Hamil- 
tonian keeps only eight bands near the Fermi level up to 15 
eV (^ 1 Ry). By the downfolding, kinetic and interaction en- 
ergies are separately renormalized into the low-energy eight 
bands and the effective Hamiltonian, where we employ the 
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Energy (eV) 

FIG. 6: (Color online) Dependence of SECI spectra (red lines) on 
scaling factor x\ (a) x = 0.8, (b) x = 1.0, and (c) x = 1.2. The closed 
blue and open blue circles represent experimental data taken from 
Ref. 31 and Ref. 32, respectively. 



CNDO model neglecting the offdiagonal part of the Coulomb 
interaction, is constructed from first principles, with the help 
of the maximally localized Wannier functions. This proce- 
dure, though several simplified treatments are employed, in 
principle, does not contain any ad hoc parameters. In the 
third stage, the Hartree-Fock method for the ground state sup- 
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plemented by the single-excitation configuration-interaction 
treatment for considering the electron-hole interactions has 
been applied to obtain electronic excitation spectra of semi- 
conductor GaAs. The spectra thus obtained have quite well 
reproduced the experimental results; the intensity and posi- 
tion for the excitonic peak are well reproduced at a quanti- 
tative level. We believe that the present model construction 
based on the downfolding offers a reliable ab initio scheme, 
where the downfolded effective Hamiltonian is capable of not 
only the ground state properties but also the excitation spectra. 

The present result opens a way of treating excitations 
such as the optical spectra by the hybrid method combin- 
ing the density functional approach and the accurate low- 
energy solver for the low-energy effective models. Beyond the 
present application to semiconductors, it would be interest- 
ing to apply this approach to excitations in strongly correlated 
electron systems such as transition metal oxides including the 
cuprates. In the present paper, we have used the Hartree-Fock 
approximation for the ground state and the single-excitation 
configuration-interaction treatment for the excitations. Opti- 
cal excitation spectra of GaAs have satisfactorily been treated 
by these approximations and the experimental results have 
been well reproduced. However, stronger electron correlation 
effects require more sophisticated low-energy solver than the 
Hartree Fock/single excitation configuration interaction treat- 
ment. For more different and challenging issues of the elec- 
tron correlation, the low-energy effective Hamiltonian may in- 
deed be treated by much more reliable low-energy solver for 
electrons in solids, such as quantum Monte Carlo methods 
for lattice Fermions,^^ path-integral renormalization group 
method,^"* and cluster extensions of the dynamical mean field 
theory.'^ 

As is well known, there are many direct ab initio 
schemes aiming at considering correlation effects; for exam- 
ple, the GW,'^ '^ transcorrelated,^' and quantum Monte Carlo 
methods.^" They are straightforward ways for approaching the 
problem compared to the present approach. However, the 
straightforward methods are faced with two serious problems: 
One is that the computational load becomes extremely heavy 
when all the electrons or even all the valence electrons are 
treated equally. The other problem is that the so far developed 
straightforward methods do not offer a sufficiently accurate 
framework if the electron correlation becomes strong such as 
in the genuine Mott insulator. The crucial point is that we need 
to treat dynamical as well as spatial correlations and fluctua- 
tions near the Fermi level in a controllable way. In the present 
stage of the computer power, such sufficient accuracies are un- 
dertaken only within simple models, which can be achieved in 
the low-energy effective model after downfolding. In fact, the 
high accuracy required from the temporal and spatial quantum 
fluctuations is important only in the low-energy region near 
the Fermi level, which justifies to restrict the high-accuracy 
treatment only in the region of low-energy excitations and thus 
only within the downfolded Hamiltonian. Within the present 
computer power, this downfolding procedure opens an avenue 
of studying highly correlated electron systems as well as ex- 
citations without relying on ad hoc parameters. By explicitly 
considering the energy hierarchy in the electronic structure. 



the first principles calculations become tractable even when 
the electron correlation is essential. 
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APPENDIX A: DERIVATION OF THE FITTING FUNCTION 
OF EQ. (10) 

Here, we describe the details of the fitting function used in 
the onsite-parameter determination. In the atomic limit, the 
onsite Hamiltonian for the decoupled site / in the home cell is 

written as 

H iJ.<f n 

= T E - 1) + Y E ^M/oiV.IO 

+ eiY^N^io, (Al) 

where N^jq and N^iq = -^^/o ^^e number operators. 
Ui and are the onsite intraorbital and interorbital Coulomb 
repulsions, respectively, e/ is a chemical potential. The on- 
site energy E(9 derived from Eq. (Al) is expressed with the 
atomic-limit wave function as 

E'c"" = {^AL\ni?\'^AL) 

= E il/J-IO - 1) + ^ E It^IOluIO 

+ e/E^c^^o, (A2) 

where we used N^^jo = (lnio \^al)- We introduce 

5q^lIo = qfiio — qio with qiQ defined as the orbital occupancy 
at the equilibrium state. At the equilibrium state, the term 
linear in Sq^io should vanish, which results in cancellation 
of the chemical potential term with the linear term from the 
Coulomb contribution. Then the quadratic energy difference 
due to the charge fluctuations is derived as 

A£™ = Y E (^9^™)' + T E -59^/0^9.70 . (A3) 
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By noting Y./.i^Qi^iof = (E^^9/^/o) - 
^ij.jiv^liJ'ioSqvio and defining the site-charge fluctua- 
tion SQio = Sq^io^ , we obtain the form 

A£;^o = Y i^Qiof + l{Ui- Ui) ^ (5g^/o<5(Z./o.(A4) 

One may specialize the charge fluctuation of one orbital be- 
cause of the crystallographycal symmetry in the system. Thus, 
the orbital index in the orbital-charge fluctuation is dropped 
and the cross term in Eq. (A4) is rewritten in terms of SQiq 
and 5qio, 

^ ^Quiohvio = 4:6qio {SQio - Sqio) ■ (A5) 

Inserting the above expression into Eq. (A4) leads to the fol- 
lowing expression for the onsite interaction energy: 



Ui 



{SQiof + 2{U'j-Ui)SQio5qio 



+ [2iUi-U'j)]i5qiof. 



(A6) 



APPENDIX B: EVALUATION OF THE MATRIX ELEMENT 

OF d/dr 

Here we describe details for the calculation of the matrix 
element {<l>rk\d / dr\<l)ak) in Eq. (32). We first rewrite it in 
terms of the maximally localized Wannier function as 



d_ 
dr 



(Bl) 



N 



(lU ij RR 

d_ 
dr 



WvjR' 



With = (}/^)T.,iCUk)Y.Re'''''\w,,R). In 

our calculation, the Wannier function Wij,i{r — R) = 
{f\'>^niR) is Stored as numerical data on the real-space grid. 



mi 



m2 

Wo 



La, 



(B2) 



where Li(= Nitti) is a cell vector defining a superlattice 
containing N{= N1N2N3) primitive cells, rui runs on the 



integer values: 0,1,- •• ,Mi — 1 with Mi being the total 
number of the grids in the ith direction (in the present case. 
Ml = M2 = M3 = 240). Since flie Wannier function satis- 
fies the foUowing periodic boundary condition. 



w^iir + Li) =w^^{r), i= 1,2,3, 



(B3) 



we can express w^i{r) in terms of the Fourier transformation 
as 



where Gl = giL\ 



[r] = ^w 

Gl 



{G 



iG, 



(B4) 



-53-^^3 withX* = {2TT/V)Lj 



Lk, V = {Li ■ L2 X X3) is the volume of the superlattice, 
and gi takes value from 1 to Mj — 1. We note that Gl is dif- 
ferent from G used in the ab initio band calculations; the for- 
mer is expressed in terms of the reciprocal-lattice vectors for 
the superlattice {L\, while the latter is written with 

the basic reciprocal lattice vectors {61, ^2, b^}. Substituting 
Eq. (B4) into Eq. (Bl) and using the translational symmetry 
for the mattix element leads to 



X Y,g^i,j{R)e'*'-^ (B5) 



R 



with 



9^ivj{R) = iV'^w*^i{GL)w^j{GL)GLe 



zGl R 



(B6) 



The actual calculation proceeds as follows: We first transform 
the real-space Wannier function w^i{r) into the reciprocal- 
space one w^i{GL) in Eq. (B4) with the algorithm of the 
fast-Fourier- transformation with radix-2, 3, and 5. Then, we 
calculate the d/dr matrix in the Wannier basis [g{R) in 
Eq. (B6)] to obtain the desired quantity [Eq. (B5)]. We note 
that this numerical procedure is a so-called Wannier interpola- 
tion scheme;"*^ we construct the Wannier functions with the ab 
initio Bloch functions in the Monkhorst-Pack /c-grid and then 
interpolate the mattix elements of d/dr at the sUghtiy shifted 
fc-grid used in the spectral calculation. 
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